Stable method and apparatus for solving S-shaped non-linear functions utilizing modified Newton-Raphson algorithms

ABSTRACT

An apparatus and method are provided for solving a non-linear S-shaped function F=f(S) which is representative of a property S in a physical system, such saturation in a reservoir simulation. A Newton iteration (T) is performed on the function f(S) at S v  to determine a next iterative value S v+1 . It is then determined whether S v+1  is located on the opposite side of the inflection point S c  from S v . If S v+1  is located on the opposite side of the inflection point from S v , then S v+1  is set to S l , a modified new estimate. The modified new estimate, S l , is preferably set to either the inflection point, S c , or to an average value between S v  and S v+1 , i.e., S l =0.5(S v +S v+1 ). The above steps are repeated until S v+1  is within the predetermined convergence criteria. Also, solution algorithms are described for two-phase and three-phase flow with gravity and capillary pressure.

TECHNICAL FIELD

The present invention relates generally to solution methods employingNewton-Raphson algorithms, and more particularly, to implicit solutionalgorithms for solving non-linear equations such as S-shaped functionsassociated with transport problems of reservoir simulators.

BACKGROUND OF THE INVENTION

Implicit solution algorithms are widely used to solve nonlinear timedependent transport problems in computational fluid dynamic (CFD)applications such as reservoir simulations. These implicit solutionalgorithms are commonly used to overcome time step size limitations ofexplicit methods. Often a solution to these transport problems for a newtime level is obtained utilizing the well-known Newton-Raphson method.

FIG. 1 schematically illustrates a simple execution of theNewton-Raphson method for solving, or at least approximating, thesolution for a function F=f(S). FIG. 1 depicts a graph of the function:F=f(S).  (1)

Solving f(S)=F=0 determines a root, i.e., where the function f(S)intersects with the horizontal axis, i.e., F=0. This intercept locationfor S is designated with the letter A. If this location, or asatisfactorily close approximation there to, can be determined, thenEqn. (1) is deemed to be solved.

The Newton-Raphson method requires making an initial guess for thesolution. The initial guess is labeled as S⁰ in FIG. 1. A correspondingpoint B on the curve F=f(S) has the coordinate (S⁰,f(S⁰)). A tangentline is drawn at point B which intercepts with the horizontal axis. Theslope of this tangent line is f′(S⁰). Algebraically, the equation of thetangent at a general point S⁰ is:F−f(S ⁰)=f′(S ⁰)(S−S ⁰)  (2)

The solution S is the horizontal coordinate of the tangent line whenF=0. Solving for this coordinate:S=S ⁰ −f(S ⁰)/f′(S ⁰)  (3)

This coordinate is then used as a second guess or approximation for thesolution of f(S)=0 and is designated by S¹. Accordingly,S ¹ =S ⁰ −f(S ⁰)/f′(S ⁰)  (4)

A better estimate of the root of Eqn. (1) can be determined by againtaking a tangent line at the coordinate (S¹,f(S¹)), i.e., point C, anddetermining its intercept with the horizontal axis. This is accomplishedusing:S ² =S ¹ −f(S ¹)/f′(S ¹)  (5)

It can be seen that S² is a better approximation than S¹ for thesolution of Eqn. (1). An even better estimate is made by taking thetangent to point D, (S²,f(S²)), which results in the intercept at S³. Ingeneral, if an initial approximation or guess for the solution of F=f(S)is S^(n) the next approximation S^(n+1) may be determined using:S ^(n+1) =S ^(n) −f(S ^(n))/f′(S ^(n))  (6)

Eqn. (6) provides the basic formula for the Newton-Raphson method. Byrepeatedly applying the Newton-Raphson Eqn. (6) to a function f(S), abetter and better approximation for the solution to Eqn. (1) can bedetermined. However, this assumes that application of the Newton-Raphsonmethod to find a solution on f(S) is unconditionally stable. This is notalways the case.

In reservoir simulation, the function to be solved is often a transportproblem which is linearized around the best estimate of the solution andsolved implicitly in each iteration. The transport problem is solved foreach cell in a reservoir model in each Newton-Raphson iteration in atypical reservoir simulator. Once a converged solution is obtained aftera number of Newton-Raphson iterations for a time step, the simulationthen proceeds with the next time step. This approach is straightforwardand unconditionally stable, if the flux function is convex (shock) asshown in FIG. 1, concave (rarefaction wave), or linear (contactdiscontinuity), which is the case for most computational fluid dynamics(CFD) problems. However, for S-shaped flux functions, which areencountered in most oil reservoir simulations, the Newton-Raphson methoddoes not converge if the time steps selected are too large. In reservoirsimulation, this problem is usually overcome by empirical time stepcontrol techniques. However, often the resulting time step size isconservative or too large, resulting in redundant computation or in timestep cuts.

Accordingly, there is a need for an implicit solution algorithm forsolving transport problems with S-shaped flux functions which isunconditionally stable, regardless of time step size. The presentinvention provides a solution to this need.

SUMMARY OF THE INVENTION

A method is provided for solving a non-linear S-shaped function F=f(S)which is representative of a property S in a physical system. TheS-shaped function has an inflection point S^(c) and a predeterminedvalue of F to which f(S) converges. An initial guess S^(v) is selectedas a potential solution to the function f(S). A Newton iteration (T) isperformed on the function f(S) at S^(v) to determine a next iterativevalue S^(v+1). It is then determined whether S^(v+1) is located on theopposite side of the inflection point S^(c) from S^(v). If S^(v+1) islocated on the opposite side of the inflection point from S^(v), thenS^(v+1) is set to S^(l), a modified new estimate. A determination isthen made as to whether S^(v+1) is within a predetermined convergencecriteria. If not, then S^(v) is set to S^(v+1) and the above steps arerepeated until S^(v+1) is within the predetermined convergence criteria.This converged value of S^(v+1) is then selected as a satisfactorysolution S to the S-shaped function F=f(S). The modified new estimate,S^(l), is preferably set to either the inflection point, S^(c), or to anaverage value between S^(v) and S^(v+1), i.e., S¹=0.5(S^(v)+S^(v+1)).Alternatively, an under-relation factor other than 0.5 may be used, suchas a value selected from 0.3 to 0.7.

Also, solution algorithms are described for two-phase and three-phaseflow with gravity and capillary pressure.

It is an object of the present invention to provide a method andapparatus utilizing an unconditionally stable modified Newton-Raphsonscheme to solve S-shaped non-linear functions to approximate solutionsfor problems based upon physical systems.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other objects, features and advantages of the presentinvention will become better understood with regard to the followingdescription, pending claims and accompanying drawings where:

FIG. 1 graphically illustrates the iterative solution of a function F=f(S) using a conventional Newton-Raphson method;

FIG. 2 is a flow diagram for one time step in a conventional reservoirsimulator wherein the operator T represents solving a linearizedtransport equation using a Newton iteration;

FIG. 3 shows a permeability field k of a 30×30 grid for a 2D test casewhich includes the locations of injector and producer wells;

FIGS. 4A, 4B, and 4C show the saturation distributions of thepermeability field k from FIG. 3 after a few initial time steps;

FIGS. 5A and 5B illustrate convergence maps of respective convex andconcave functions f(S);

FIGS. 6A and 6B depict S-shaped flux functions and FIGS. 6C and 6Dillustrate corresponding convergence maps made utilizing a conventionalNewton-Raphson scheme;

FIGS. 7A and 7B illustrate convergence maps obtained using a firstmodified Newton-Raphson scheme, made in accordance with the presentinvention, for solving the S-shaped functions shown in FIGS. 6A and 6B;

FIG. 8 is a flow diagram of a second modified Newton-Raphson scheme forone time step wherein solving the linearized transport equation isrepresented by operator T;

FIGS. 9A and 9B display convergence maps using a second embodiment of amodified Newton-Raphson scheme for solving the S-shaped functions shownin FIGS. 6A and 6B;

FIGS. 10A, 10B, 10C, 10D, 10E, 10F, 10G, 10H, and 10I, respectively,show saturation fields after time 0.5 Pore Volume Injected (PVI) forthree ratios of viscosities μ₁/μ₂ and three time step sizes;

FIG. 11 illustrates convergence histories after time 0.5 PVI and μ₁/μ₂=1using first and second modified Newton schemes and using a Newton schemewith constant under-relaxation;

FIG. 12 is a flow diagram with a non-constant velocity field for onetime step wherein solving a pressure equation is represented by theoperator P and solving the linearized transport equation is representedby the operator T;

FIG. 13 shows the convergence histories of one time step for threesimulations, two performed with modified Newton-Raphson schemes and athird using constant under-relaxation during all iterative steps;

FIG. 14 is a flow diagram of steps taken in a preferred embodiment ofthe present invention to solve an S-shaped function, F=f(S), for aphysical property S, such as saturation;

FIG. 15 graphically illustrates the iterative solution of an S-shapedfunction F=f(S) using a modified Newton-Raphson scheme in accordancewith the present invention;

FIGS. 16A and B illustrate fractional flow curves for phase 1:

(A) m=1 and (B) m=0.2;

FIGS. 17A and B show fractional flow curves for phase 2: (A) m=1 and (B)m=0.2;

FIGS. 18A and B depict fractional flow between cells i and n: in theupper figure, cell i is in upwind and in the lower figure cell n is inupwind, (A) no counter-current flow; (B) counter-current flow in cell i(C) counter-current flow in cell n; and (D) counter-current flow in bothcells i and n;

FIGS. 19A and B show oil phase (lighter fluid) fractional flow curveswith counter-current and co-current flows;

FIGS. 20A and B illustrate fractional flow surfaces with no gravityeffect for (A) phase 1 and (B) phase 2: m₁=0.2, m₃=10, andC_(g1)=C_(g2)=0;

FIGS. 21A and B show fractional flow surfaces with gravity effect for(A) phase 1 and (B) phase 2: m₁=0.2, m₃=10; C_(g1)=5, C_(g2)=10;

FIGS. 22A and B depict fractional flow surfaces with no gravity effects(A) phase 1 and (B) phase 2: m₁=0.2, m₃=10. C_(g1)=0, C_(g2)=0;

FIGS. 23A and B depict fractional flow surfaces with small gravityeffects (A) phase 1 and (B) phase 2: m₁=0.2, m₃=10. C_(g1)=0.5 andC_(g2)=1.0;

FIGS. 24A and B depict fractional flow surfaces with large gravityeffects (A) phase 1 and (B) phase 2: m₁=0.2, m₃=10. C_(g1)=5 andC_(g2)=10;

DETAILED DESCRIPTION OF THE INVENTION

I. Theoretical Background and Numerical Studies for Reservoir Simulation

Consider the following hyperbolic equation: $\begin{matrix}{{\frac{\partial S}{\partial t} + {\nabla({fu})}} = {0\quad{on}\quad\Omega}} & (7)\end{matrix}$

Eqn. (7) is a model equation which represents the transport of one phasein a two phase subsurface flow problem. The dependent variable S is thesaturation of that phase (0≦S≦1), f is a function of S with 0≦f≦1, and uis a divergence free total velocity field. For simplicity, but withoutloss of generality, u is considered constant, and it is assumed that∂f/∂S≧0.

The nonlinear transport problem of Eqn. (7) may be solved iterativelyusing a conventional Newton-Raphson method as outlined in the flowdiagram of FIG. 2. This procedure results in the solution S^(n+1) at anew time level. An initial guess S^(v) for the solution to an S-shapedfunction is made. This guess may be the saturation value S^(n) from theprevious time step. In each iteration v, the linearized transportequation is preferably solved for each cell in a reservoir model orsubmodel (represented by the operator T) which results in an updatedsaturation field S^(v+1). Once the maximum absolute saturation changefrom one iteration to the next does not exceed an appropriately definedthreshold ε, this nonlinear loop is considered converged and S^(v+1)constitutes the solution S^(n+1) at the new time level.

A 2D discretization of Eqn. (7) is studied: $\begin{matrix}{{{S\overset{v + 1}{i,j}} - {S\overset{n}{i}}},{{j + {\underset{\underset{c_{x}}{︸}}{u_{x}\frac{\Delta\quad t}{\Delta\quad x}}\left( {\overset{n + 1}{f_{i,j} -}\overset{n + 1}{f_{{i - 1},j}}} \right)} + {\underset{\underset{c_{y}}{︸}}{u_{y}\frac{\Delta\quad t}{\Delta\quad y}}\left( {\overset{n + 1}{f_{i,j} -}\overset{n + 1}{f_{i,{j - 1}}}} \right)}} = 0}} & (8)\end{matrix}$which is based on implicit Euler and first order upwinding. Forsimplicity, Eqn. (8) is written for the case in which both components ofthe velocity vector u={u_(x) u_(y)}^(T) are greater than or equal tozero and an orthogonal, uniformly spaced grid is assumed. Each grid cellcan be identified by an index pair i, j and the grid spacing in the twocoordinate directions is defined by Δx and Δy. The superscripts n and vdenote the old time and iteration levels, respectively, and Δt is timestep size. In order to solve Eqn. (8), the flux f^(n+1) in Eqn. (8) isapproximated by the linearization $\begin{matrix}{{f^{n + 1} \approx f^{v + 1}} = {f^{v} + {\underset{\underset{f^{\prime\quad v}}{︸}}{\left( \frac{\partial f}{\partial S} \right)^{v}}\left( {S^{v + 1} - S^{v}} \right)}}} & (9)\end{matrix}$

Finally, an approximation S^(v+1) of S^(n+1) is obtained by solving thelinearized system:S _(i,j) ^(v+1) +c _(x)(f′ _(i,j) ^(v) S _(i,j) ^(v+1) −f′ _(i−1,j) ^(v)S _(i−1,j) ^(v+1))+c _(y)(f′ _(i,j) ^(v) S _(i,j) ^(v+1) −f′ _(i,j−1)^(v) S _(i,j−1) ^(v+1))=S _(i,j) ^(n) −c _(x)(f _(i,j) ^(v) −f _(i−1,j)^(v))−c _(y)(f _(i,j) ^(v) −f _(i,j−1) ^(v))+c _(x)(f′ _(i,j) ^(v) S_(i,j) ^(v) −f′ _(i−1,j) ^(v) S _(i−1,j) ^(v))+c _(y)(f′ _(i,j) ^(v) S_(i,j) ^(v) −f′ _(i,j−1) ^(v) S _(i,j−1) ^(v))  (10)which is obtained by substituting Eqn. (9) into Eqn. (8). It isimportant to enforce 0≦S^(v+1)≦1 every time Eqn. (10) is solved.

For numerical studies presented in this specification, a 30×30 grid ofcells is used for modeling. FIG. 3 shows a permeability field k andinjector and producer well locations employed to compute the totalvelocity u=−λΔp by solving the elliptic equation:∇·(λ∇p)=−qonΩ  (11)for the pressure p. Here, λ is the total mobility$\lambda = \left( {\frac{k_{r\quad 1}}{\mu_{1}} + \frac{k_{\quad{r\quad 2}}}{\mu_{2}}} \right)$where k_(rj) is the relative permeability of phase j and μ_(j) is theviscosity of the phase. The source term q is 1 at the producer well and−1 at the injector well and S≡0 is the initial condition for thesaturation. For the case of a simple quadratic function, an exemplaryS-shaped fractional flow function can be represented by $\begin{matrix}{{f(S)} = \frac{\frac{S^{2}}{\mu_{1}}}{\frac{S^{2}}{\mu_{1}} + \frac{\left( {1 - S} \right)^{2}}{\mu_{2}}}} & (12)\end{matrix}$where μ₁ and μ₂ are the viscosities of the two phases (the injectedphase is labeled with the subscript 1). In order to demonstrate that theimplicit solution algorithm explained above, with respect to FIG. 2, isunstable for time steps which are too large, results are presented forthree numerical experiments for which μ₁/μ₂ is chosen to be 1, 10, and0.1, respectively. FIG. 4 shows the saturation distributions after a fewinitial time steps. Although the time steps are relatively small interms of pore volume injected (PVI), the results show strong oscillatorybehavior characterized by the checkerboard-like pattern of thesaturation front. The nonlinear loop of the solution algorithm of FIG. 2does not converge and the number of iterations is limited to 200. Theorigin of this instability, or oscillatory behavior, will now beinvestigated and a modification will be described which results in anunconditionally stable scheme for use of the Newton-Raphson method.

A key to understanding the stability issues of nonlinear transportproblems lies in the linearization of Eqn. (9) of the flux function.First, the nonlinearity of the transport problem f(S) is isolated.Second, a determination is made of those starting values, in the range0≦S⁰≦1, for which an iterative scheme fulfills the convergenceconditions: $\begin{matrix}{{{\lim\limits_{v->\infty}{f\left( S^{v} \right)}} = F},} & (13) \\{{{f\left( S^{v} \right)} + {{f^{\prime}\left( S^{v} \right)}\left( {S^{v + 1} - S^{v}} \right)}} = {F.}} & (14)\end{matrix}$of f(S) around S^(v). The Newton iteration equation is obtained asfollows: $\begin{matrix}{S^{v + 1} = {S^{v} + {\frac{F - {f\left( S^{v} \right)}}{f^{\prime}\left( S^{v} \right)}.}}} & (15)\end{matrix}$

In addition, 0≦S≦1 is enforced after each iteration. Now thecombinations of S⁰ and F for which condition (13) is met will beinvestigated. The condition (13) is always met if the function f is alinear, convex, or concave function of S. Since the linear case istrivial, only examples for the convex and concave cases will be shownwhich are, respectively, represented by f(S)=S² and f(S)=1−(1−S)².

FIGS. 5A and 5B show convergence maps of these convex and concavefunctions. The convergence maps are created by selecting an initialvalue of S⁰ and a target value F for the flux function f. Then thenumber of iterations is determined which are required to converge withina predetermined tolerance to target value F. The horizontal axisrepresents the start value S⁰ and the vertical axis the target value F.This process is repeated for a large number of start values 0≦S⁰≦1 attarget value F. The process is then repeated for a number of targetvalues of 0≦F≦1. This determination of the number of iterations requiredto achieve convergence for each value (S, F) is then plotted to producea convergence map.

Fast convergence is represented by darker shades and white color showsthe regions in the parameter space for which the iteration Eqn. (15)does not converge. There is no white region present in FIGS. 5A and 5Band therefore the iteration scheme is unconditionally stable for convexand concave functions f(S).

The situation is different, however, if the flux function f is anS-shaped function. An S-shaped function is characterized by having aconvex portion and a concave portion which are joined at an inflectionpoint S^(c) as seen in FIGS. 6A and 6B. The inflection point S^(c) ischaracterized as being the point about which the second derivatives ofthe function f(S), f″(S), changes sign.

The upper plots in FIGS. 6A and 6B show f defined by Eqn. (12) forμ₁/μ₂=1 and μ₁/μ₂=0.1, respectively. The corresponding convergence mapsare shown in FIGS. 6C and 6D. The existence of large, white regionsclearly shows that the iteration scheme applied to the flux function fof Eqn. (12) is not unconditionally stable for these values of μ₁/μ₂.Convergence can be achieved if the initial guess for S⁰ is close to thefinal solution, i.e., if |f(S⁰)−F| is not too large. This is consistentwith limiting the time step size or the maximum saturation change, if atransport problem is solved.

Convergence maps in FIGS. 6C and 6D show that the iteration Eqn. (15)converges for all possible target values F, if S⁰=S^(c), where S^(c) isthe inflection point of f. This finding suggests setting S^(v+1)=S^(c),if the second derivatives of f(S^(v)) and f(S^(v+1)) have oppositesigns. If this is done at the end of each Newton iteration, a modifiedNewton-Raphson scheme is produced which is unconditionally stable forS-shaped flux functions.

FIGS. 7A and 7B show convergence maps for the S-shaped flux functions ofFIGS. 6A and 6B. However, this time a first modified Newton-Raphsonscheme is applied of setting S^(v+1)=S^(c) if the second derivatives off(S^(v)) and f(S^(v+1)) have opposite signs. This can be verified bycalculating:f″(S ^(v))f″(S ^(v+1))<0  (16)

If the product of second derivatives is negative, then S^(v) and S^(v+1)are on opposite sides of the inflection point S^(c).

The horizontal axis represents the start value S⁰ and the vertical axisis the target value F. Dark regions indicate fast convergence, and noconvergence is achieved for the parameter pairs in white regions. Notethat no white regions are present in FIGS. 7A and 7B and thus this firstmodified Newton-Raphson scheme is unconditionally stable. Theapplication of this modification to the solution algorithm for thetransport problem is straightforward and allows for arbitrarily largetime steps.

It can be difficult to determine the exact location of the inflectionpoint S^(c). This is particularly true if a system of transportequations is solved, where the fractional flow functions depend onmultiple variables. In some instances, look-up tables may be used torepresent the fraction flow curve f(S) and the inflection points are notreadily computed analytically. Also, finding the inflection point S^(c)of the complex fraction flow functions can be computationally expensive.

Therefore, an alternative, and more preferred, second modified scheme ispresented which is easier to implement in a realistic setting. Thissecond scheme differs from the first scheme in the way S^(v+1) ismodified at the end of each Newton iteration where S^(v) and S^(v+1) arelocated on opposite sides of the inflection point S^(c) of curve f.Instead of setting S^(v+1) to S^(c), a selective under-relaxation methodis applied. For example, S^(v+1) may be replaced by 0.5(S^(v)+S^(v+1)),if the sign of the second derivative, i.e., f″(S)=∂²f/∂S², is not thesame at the old and new iteration levels.

Other levels of relaxation may also be used. That is, theunder-relaxation can be generalized by:S ^(l) =αS ^(v)+(1−α)S ^(v+1)  (17)where S^(l)=a new estimate for S; and α=an under-relaxation coefficient.

A flow diagram of this second enhanced scheme is depicted in FIG. 8.Since it is relatively easy to determine the second derivatives of f,this second alternative scheme can easily be implemented in an existingoil reservoir simulator as compared to the first modified scheme whereS^(v+1) is set to S^(c) when S^(v) and S^(v+1) are on opposite sides ofS^(c).

Analogue to the convergence maps shown in FIGS. 6C and 6D and FIGS. 7Aand 7B, the corresponding convergence maps for the second alternativemodification using selective under-relaxation are shown in FIGS. 9A and9B. The horizontal axis represents the start value S⁰ and the verticalaxis represents the target value F. Dark regions indicate fastconvergence, and no convergence is achieved for parameter pairs in whiteregions. As with the convergence maps of in FIGS. 7A and 7B, there areno white regions showing that this second alternative modification leadsto an unconditionally stable scheme.

Both the first and second modifications work as predicted, if they areapplied to the solution algorithm for nonlinear transport equations. Thefirst and second modified schemes are almost equally efficient.Moreover, if the time step size is within the stable range for theunmodified or conventional Newton-Raphson scheme, no compromise inefficiency due to the modifications of this invention could be observedin numerical studies. On the other hand, efficiency is heavilycompromised, if stability is achieved by applying under-relaxationeverywhere at the end of each Newton iteration rather than selectivelyas in the modified schemes described above. Numerical results arepresented in FIGS. 10A-I for the test case described above usingS-shaped fractional flow functions defined by Eqn. (12). All simulationswere performed with the first and second modified schemes and it isdemonstrated that arbitrary large time steps can be applied. FIGS.10A-I, respectively, show saturation fields obtained after time=0.5 PVIwith μ₁/μ₂=1 (top row—FIGS. 10A, 10B, and 10C), μ₁/μ₂=10 (middlerow—FIGS. 10D, 10E, and 10F), and μ₁/μ₂=0.1 (bottom row—FIGS. 10G, 10H,and 10I). The results in the left column of these plots were computedusing 100 small time steps, or Δt=0.005 PVI, for which an unmodifiedNewton-Raphson scheme proved to be stable as well. Note that these timesteps are half the size of those used to obtain the non-convergingresults shown in FIG. 4. The center and right columns of plots showresults computed with Δt=0.1 PVI and Δt=0.5 PVI, respectively.

FIG. 11 shows the results of using the above described first and secondmodified Newton-Raphson schemes wherein S^(v+1) is selectively set toS^(c) and 0.5(S^(v)+S^(v+1)), respectively, only when the signs of thesecond derivatives, f″(S^(v)) and f″(S^(v+1)), change in an iterationstep. Further, a third scheme, was investigated wherein under-relaxationwas applied at the end of every iteration, i.e., a naïve under-relationscheme. The vertical axis shows maximum saturation change in aniterative step and the horizontal axis shows number of iterations. FIG.11 is plotted for the case where the time step size was 0.5 PVI and theratio of viscosities μ₁/μ₂ was one. For all simulations performed, boththe first and second modified schemes were tested and no noticeabledifference in terms of required iterations could be observed. While theconvergence rates of the first two schemes are essentially the same, thethird scheme with naïve under-relaxation is shown to be much lessefficient.

The modified Newton-Raphson schemes can also be applied as a buildingblock for a solution algorithm, if mobility factor λ is a function ofthe saturation. This is of particular interest for the oil industry. Forthe quadratic relative permeability case, mobility factor λ is definedas: $\begin{matrix}{\lambda = {k\left( {\frac{S^{2}}{\mu_{1}} + \frac{\left( {1 - S} \right)^{2}}{\mu_{2}}} \right)}} & (18)\end{matrix}$

-   -   where k is permeability;    -   S=saturation of a first phase;    -   μ₁=viscosity of the first phase; and    -   μ₂=viscosity of a second phase.

As a consequence of λ being a function of the saturation, the velocityfield u is no longer constant and has to be updated by solving thepressure Eqn. (18), if λ changes. This algorithm is outlined by the flowdiagram in FIG. 12. Note that the flow diagram of FIG. 8 reappears inthe dashed box within the loop required to update the total velocityfield u. This converged solution is identical to the solution obtainedby a fully implicit scheme, where the pressure and transport equationsare solved simultaneously.

FIG. 13 shows the convergence histories of one time step for twosimulations, one performed with the modified Newton-Raphson solutionalgorithm using selective under-relaxation and the other one with amethod using simple or naïve under-relaxation. The test case was thesame as for the previous studies, the time step size was 0.5 PVI and theviscosity ratio μ₁/μ₂=1. With the modified algorithm, the linearizedtransport equation has to be solved approximately 60 times and theelliptic pressure equation only nine times (number of peaks in theconvergence plot). Simple or naïve under-relaxation on the other handproved to be much less efficient.

Exemplary Steps for Solving an S-Shaped Function for a Physical Property

The present invention includes a method and apparatus for stably solvingnon-linear S-shaped functions, F=f(S), to determine a property Srepresentative of a physical system. FIG. 14 illustrates a preferredexemplary flow diagram of steps which may be used to carry out theinvention in accordance with principles described in the section above.FIG. 15 graphically illustrates a modified Newton scheme.

In this modified Newton-Raphson algorithm, the value for a nextiteration, S^(v+1), is set to a new estimate or stability limit S^(l) inthe event the signs of the second derivatives, f″(S^(v)) and f″(S^(v+1))of the S-shaped function change between successive iterations. Themodified new estimate, or stability limit, S^(l) is ideally selectedsuch that the modified Newton-Raphson methodology is guaranteed tostably converge.

A first step 110 is to obtain an S-shaped function f(S) which isrepresentative of a property S in a physical system. For example, aphysical system of particular interest is that of fluid flow in asubterranean reservoir. The S-shaped function f(S) may be expressed interms of a mathematical function, such as Eqn. (12) above, or else byway of a look-up table to capture the essence of the S-shaped function.As is well known to those skilled in the art of reservoir simulation, anS-shaped function, such as a fractional flow function, may be derivedfrom laboratory tests based on relative permeability measurements.

Those skilled in the art of solving non-linear problems for physicalsystems will appreciate that the present invention, utilizing modifiedNewton-Raphson schemes, may be used to solve S-shaped functions whichcharacterize properties of physical systems other than reservoirsimulation.

The S-shaped function f(S) is characterized by having at least a convexportion and a concave portion joined at an inflection point S^(c). Thefunction f(S) is to be solved to determine a physical property S whichsatisfies the condition F=f(S), where F is a target value. For example,in a reservoir simulation, when solving an S-shaped fractional flowfunction f(S) for a saturation value S in a grid cell in a linearizedtransport problem, the flow flux value F from a previous time level isideally used as the target value F. Alternatively, target value F may beestablished by non-linear iteration in numerical simulation.

An initial value or guess S^(v) is selected in step 120 as a possiblesolution for the physical property S which satisfies the conditionF=f(S). This initial value for S^(v) is preferably established by aninitial condition or the results from a preceding time step. Althoughless preferred, the initial value for S^(v+1) may also be made by aninterpolation method or randomized method.

A Newton iteration, as suggested by the operator (T) in FIG. 8, isperformed in step 130 on the function f(S) utilizing the value S^(v) ofthe present iteration to determine a next iteration value S^(v+1). ThisNewton iteration corresponds to applying Eqn. (15) to the S-shapedfunction f(S) to determine the next iterative value S^(v+1). In FIG. 15,this initial guess is S⁰.

An optional check may be made in step 140 to see whether S^(v+1) iswithin an acceptable bound. For example, if the denominator in anS-shaped function, such as a function similar to that of Eqn. (12), werenear zero, the estimated value for S^(v+1) could be quite large. IfS^(v+1) exceeds certain boundary limits, such as saturation S being lessthan zero or greater than one, then S^(v+1) may be set in step 150 tovalues such as 0.0 or 1.0 to maintain the physical property S, withinrealistic bounds. In the example of FIG. 15, the intercept of thetangent to point (S⁰, f(S⁰)) with the horizontal line representing flowflux value F is seen not to occur within the upper saturation limit,S^(b)=1.0. Accordingly, S^(v+1) in the case of the first iteration inFIG. 15, is set to the upper bound S^(b)=1.0.

A determination is then made in step 160 as to whether S^(v) and S^(v+1)reside on opposite sides of the inflection point S^(c) of the S-shapedfunction f(S). The preferred manner of making such a determination is tocheck to see whether the signs of the second derivative of f(S), i.e.,f″(S), at S^(v) and at S^(v+1) are opposite one another. This can beverified utilizing Eqn. (16) to calculate the product of the f″(S^(v))and f′(S^(v+1)).

If the product of second derivatives is negative, then S^(v) and S^(v+1)are on opposite side of the inflection point. Another, although lesspreferred, manner of making this determination is to actually calculatethe location of the inflection point S^(c) and then see whether S^(v)and S^(v+1) are on the same side of S^(c). The inflection point S^(c)may be calculated by solving f″(S)=0. This may be done analytically ifthe function f(S) is sufficiently simply or else numerically where f(S)is complex. In the case of using a look-up table to define the S-shapedfunction f(S), an initial table of values for f(S), f′(S) and f″(S) maybe generated. The look-up tables may then be examined to determine thelocation of the inflection point S^(c). In FIG. 15, values of S⁰ andS^(b) are located on opposite sides of the inflection point S^(c).

If S^(v) and S^(v+1) are on opposite sides of the inflection pointS^(c), then to guarantee unconditional stability in solving for S,S^(v+1) is set to S^(l), a new estimate or stability limit, in step 170.This modified estimate S^(l) is selected to guarantee convergence of theNewton-Raphson iterations. In a first preferred embodiment, the modifiedestimate S^(l) is set equal to the inflection point S^(c). In a secondand more preferred manner as shown in FIG. 15, an under-relaxationmethod is used to calculate S^(l). Ideally, S^(l) is calculated inaccordance with the following equation:S ^(l) =αS ^(v)+(1−α)S ^((v+1))  (19)where 0<α<1.

Preferably the range of α is 0.3<α<0.7, even more preferably 0.4<α<0.6,and most preferably, α=0.5. With α=0.5, S^(v+1) or S^(l) is set to0.5(S^(v)+S^(v+1)).

If it is determined that S^(v) and S^(v+1) are located on the same sideof inflection point S^(c), then S^(v+1) is preferably left unalteredrather than being set to the modified estimate S^(l), i.e., S^(c) orS^(l)=0.5(S^(v)+S^(v+1)).

In step 180, a determination is made whether the current iterative valueS^(v+1) meets a predetermined tolerance condition. For example, it maybe required that the difference between consecutive iterative values ofS^(v) and S^(v+1) be sufficiently small. Mathematically, this may beexpressed as:ε>|S ^(v+1) −S ^(v)|  (20)where ε is a predetermined tolerance limit.

Alternatively, it may be required that:ε>|F−f(S ^(v+1))|  (21)in order to consider the solution S^(v+1) to be sufficiently close tothe exact solution.

If S^(v+1) is not considered a satisfactory solution, then S^(v) is setto S^(v+1) in step 190 and steps 130-180 are continued until a value ofS^(v+1) converges to within the predetermined tolerance limit. FIG. 15demonstrates how successive estimates for saturation S¹, S² and S³,etc., converge toward the exact solution for saturation S. Whensufficient convergence is reached, the physical value S is considered asbeing a satisfactory solution. In the case of a reservoir simulator, thesaturation S in grid cells for the next time level S^(n+1) can then beset to S^(v+1) after convergence is reached for each of the grid cellsover which the modified Newton iterations are being performed.

The present invention further includes an apparatus or program storagedevice readable by a machine, tangibly embodying a program ofinstructions. These instructions are executable by the machine toperform method steps for stably solving a non-linear S-shaped function,F=f(S), which is representative of a property S in a physical system.The S-shaped function has an inflection point S^(c) and a predeterminedvalue of F to which f(S) converges. The method comprises the followingsteps. An initial value S^(v) for the possible solution to F=f(S) isselected. A Newton iteration is performed on the function f(S) at S^(v)to determine a next iterative value S^(v+1). A determination is thenmade as to whether S^(v+1) is located on the opposite side of theinflection point S^(c) from S^(v). If S^(v) and S^(v+1) are on oppositesides, then S^(v+1) is set equal to S^(l), a modified new estimate. Nexta check is made to determine whether S^(v+1) is a satisfactory solutionfor the property S by determining if S^(v+1) is within a predeterminedconvergence criteria. If S^(v+1) is not satisfactorily converged, thenS^(v) is set equal to S^(v+1). The above steps are repeated untilS^(v+1) is within the predetermined convergence criteria. Whenconvergence is reached, S^(v+1) is selected as a satisfactory solutionto the problem F=f(S).

II. Algorithm for Two-Phase Flow with Gravity and Capillary Pressure

The effects of gravity and capillary pressure on the numerical stabilityof the Newton-Raphson method will be examined for the non-lineartransport equation for two-phase flow. Gravity and/or capillary pressureeffects in reservoir simulation can incur counter current flow in cellsof a reservoir. Such counter current flow can pose severe numericalinstability in conventional reservoir simulators. By extending themethod described above in section I, an unconditionally stable algorithmwill be derived for the non-linear transport equation for two-phase flowwith gravity and capillary pressure.

First, the gravity effect will be investigated on fractional flowcurves. Then, fractional flow curves will be characterized by acharacteristic parameter. An unconditionally stable algorithm fortwo-phase flow with gravity is then described. A modification is alsoprovided to the unconditionally stable algorithm to include the effectof capillary pressure.

Gravity Effect

Darcy's law for two phase flow with gravity can be written as$\begin{matrix}{u_{1} = {\frac{{kk}_{r\quad 1}}{\mu_{1}}\left( {{\nabla p} - {\rho_{1}g{\nabla D}}} \right)}} & (22) \\{u_{2} = {\frac{{kk}_{r\quad 2}}{\mu_{2}}\left( {{\nabla p} - {\rho_{2}g{\nabla D}}} \right)}} & (23)\end{matrix}$

Here, k is permeability, u is phase velocity, k, is relativepermeability, μ is viscosity, D denotes the depth and ρ is density. Thesubscripts denote phase properties.

From Eqns. (22) and (23), it can be readily derived that:u ₁ =f ₁(S)u−f ₃(S)C _(g)  (24)u ₂ =f ₂(S)u+f ₃(S)C _(g)  (25)where S denotes the saturation of phase 1 (i.e., water). The totalvelocity is defined as $\begin{matrix}{u = {u_{1} + u_{2}}} & (26) \\{f_{1} = \frac{k_{r\quad 1}}{k_{r\quad 1} + {mk}_{r\quad 2}}} & (27) \\{f_{2} = \frac{{mk}_{r\quad 2}}{k_{r\quad 1} + {mk}_{r\quad 2}}} & (28) \\{f_{3} = \frac{{mk}_{r\quad 1}k_{r\quad 2}}{k_{r\quad 1} + {mk}_{r\quad 2}}} & (29)\end{matrix}$

Here, m is the viscosity ratio, μ₁/μ₂.

The velocity due to gravity force is characterized by variable, C_(g),$\begin{matrix}{C_{g} = {\frac{{kg}\left( {\rho_{1} - \rho_{2}} \right)}{\mu_{1}u_{c}}{\nabla D}}} & (30)\end{matrix}$where u_(c) is a characteristic velocity and D is the depth. Note thatvelocities u, u₁ and u₂ in Eqns. (24) and (25) are alsonondimensionalized with respect to u_(c).

It is instructive to examine fractional flow curves from Eqns. (24) and(25). In FIGS. 16 and 17 the fractional flows are plotted for phases 1and 2, respectively, $\begin{matrix}{F_{1} = \frac{u_{1}}{u}} & (31) \\{{F_{2} = \frac{u_{2}}{u}},} & (32)\end{matrix}$with quadratic relative permeabilities,k _(r1)=S₁ ²  (33)k _(r2) =S ₂ ²=(1−S ₁)².  (34)

The fractional flow curves are monotonic if 0<mC_(g)≦1, whereas it canbecome a multi-valued function for mC_(g)>1 or mC_(g)<0. Note that thefractional flow can be smaller than zero or larger than one, dependingon the directions of gravity force and velocity. This can causenumerical difficulty if the Newton nonlinear iteration does not take theactual pattern of theses curves into account. Note that the fractionalflow curve from negative C_(g) and positive u is equivalent to that frompositive C_(g) and negative u. In practical simulation with fixedcoordinates, C_(g) can be kept always positive, but the mean velocitycan be either positive or negative. Note that without loss ofgenerality, the case with negative mean velocity is equivalent to thecase of positive mean velocity with negative C_(g) by selecting anappropriate coordinate to make the mean velocity positive.

If oil (lighter fluid) equation is chosen as a primary equation, fromEqn. (25) the saturation where the fractional flow of oil phase becomeszero S*₀ or one S*₁:S* ₀=1−√{square root over (|u|/C _(g))} for u<0, and  (35)S* ₁=1−√{square root over (|u|/mC _(g))} for u>0.  (36)

If mC_(g)/|u|≦1 for u>0, the fractional curve becomes a monotonic curve;whereas if u<0 or mC_(g)/|u|>1 for u>0, the fractional flow curveincludes a domain where the fractional curve becomes multi-valued (nonmonotonic). When the fractional flow curves has a minimum or maximumpoint, then there exists a reflection point which the saturation doesnot cross over in the non-linear iteration (see FIG. 19). Thisreflection point can be calculated analytically if the relativepermeabilities are quadratic as in Eqns. (33) and (34). $\begin{matrix}{0 = {{mS}_{r}^{*3} - \left( {1 - S_{r}^{*}} \right)^{3} + {\frac{u}{C_{g}}.}}} & (37)\end{matrix}$

The cubic equation can be solved analytically and always yields at leastone real root: $\begin{matrix}{{S_{r}^{*} = {\alpha + \left( {{- \frac{q}{2}} + \frac{\sqrt{q^{2} + {4p^{3}}}}{2}} \right)^{1/3} + \left( {{- \frac{q}{2}} - \frac{\sqrt{q^{2} + {4p^{3}}}}{2}} \right)^{1/3}}};} & (38) \\{{\alpha = \frac{1}{1 + m}};} & (39) \\{{p = {\frac{1}{3}\left( {{3\alpha^{2}} + {2b\quad\alpha} + c} \right)}};} & (40) \\{{q = \left( {\alpha^{3} + {b\alpha}^{2} + {c\quad\alpha} + d} \right)};} & (41) \\{{b = {- \frac{3}{1 + m}}};} & (42) \\{{{c = \frac{3}{1 + m}};}{and}} & (43) \\{d = {- {\frac{1 - {u/C_{g}}}{1 + m}.}}} & (44)\end{matrix}$

where

The saturation can be divided into two domains: (1) S>S*_(r) andS<S*_(r). In general the saturation does not cross over the reflectionpoint, but if the cell is next to the boundary with no-flow or next tothe cell with the saturation of the other domain, the saturation shouldbe allowed to cross over. The fractional flow curve does not apply dueto the presence of the boundary that allows the fluid to accumulatewithout counter-current flow.

Finite Difference Approximation

The fractional flow curves in FIGS. 16 and 17 indicate that two phaseflow can be co-current and counter-current, depending on thecharacteristic variable Cg and the viscosity ratio, m, and saturation S.As seen in FIG. 18, when the fractional flow in cell i and theneighboring cell n does not incur counter-current flow, a generalup-wind finite difference approximation can be applied withoutambiguity. However, if counter-current flow occurs in either cell i orcell n, as shown in FIG. 18, the upstream is phase dependent.

The phase split from the total velocity calculated in the pressuresolution, as a result, entails a detailed information of co-current orcounter-current flow between cells i and n. In FIG. 18 are depicted allthe possible cases that can occur between cell i and its neighboringcell n. It is assumed that the total velocity from the pressuresolution, u_(f), is the velocity at the interface of cells i and n andthe total velocities inside cells i and n are identical, say, u.

In cell i, the phase velocities that are co-current and counter-currentare denoted with the total velocity u by u_(p) ^(i) and u_(r) ^(i),respectively. These phase velocities can be obtained from Eqn. (24) or(25). The equations can be expressed asu _(p) ^(i) =f(u,S ^(i) ,C _(g) ,m),  (45)u _(r) ^(i) =u−u _(p) ^(i).  (46)

Similarly, the phase velocities in cell n are obtained byu _(p) ^(n) =f(u,S ^(n) ,C _(g) ,m),  (47)u _(r) ^(n) =u−u _(p) ^(n).  (48)

If cell i is upstream,u _(f) =u _(p) ^(i) +u _(r) ^(n)  (49)and if cell n is in upstream,u _(f) =u _(p) ^(n) +u _(r) ^(i).  (50)

Note that it is assumed that a particular phase velocity can beco-current or counter-current depending on the flow and gravitydirections. Furthermore, counter-current flow does not occur if0<mC_(g)≦1.

In the sequential implicit algorithm, the saturation equations areiteratively solved with the total velocity fixed to preserve theconservation of mass (i.e., divergence free). This algorithm can bestraightforwardly implemented by the total velocity and fractional flowformulation. However, for the case of counter-current flow, the originalfractional flow that is dependent on a single saturation has to bemodified in order to maintain mass conservation. If the phase velocityu; is determined from the upstream α_(i) cell, the fractional flow ofthe cell α_(i) can be expressed as: $\begin{matrix}{u_{i}^{\alpha_{i}} = {k\frac{k_{ri}^{\alpha\quad i}}{\mu_{i}^{\alpha\quad i}}\left( {{\Delta p} - {\rho_{i}g{\nabla D}}} \right)}} & (51)\end{matrix}$

Then the total velocity can be written as: $\begin{matrix}{u = {{\sum\limits_{l}u_{l}^{\alpha\quad l}} = {{k\quad\Delta\quad p{\sum\limits_{l}\frac{k_{r\quad l}^{\alpha_{l}}}{\mu_{l}^{\alpha_{l}}}}} - {{kg}{\nabla D}{\sum\limits_{l}{\frac{k_{rl}^{\alpha\quad l}\rho_{l}g}{\mu_{l}^{\alpha\quad l}}\left( {{\Delta\quad p} - {\nabla D}} \right)}}}}}} & (52) \\{u_{i}^{\alpha_{i}} = {{\frac{\lambda_{i}^{\alpha_{i}}}{\lambda_{T}}u} + {{kg}{\nabla D}\frac{\lambda_{i}^{\alpha_{i}}}{\lambda_{T}}{\sum\limits_{l}{\lambda^{\alpha_{l}}\left( {\rho_{l} - \rho_{i}} \right)}}}}} & (53) \\{\lambda_{l}^{\alpha_{l}} = \frac{k_{r\quad l}^{a_{l}}}{\mu_{l}^{\alpha_{l}}}} & (54) \\{\lambda_{T} = {\sum\limits_{l}{\lambda_{l}^{\alpha_{l}}.}}} & (55)\end{matrix}$Stability Check

The stability issues of nonlinear transport problems lie in thelinearization of the flux function. Note that the saturation equation isnow a function of the cell saturation and all the upstream-cellsaturations. The numerical schemes have been examined as to whether theschemes can converge from the initial saturation estimate, 0≦S⁰≦1, forall cells, to the final saturation S^(v) that satisfies the convergencecriterion: $\begin{matrix}{{{\lim\limits_{v->\infty}{F\left( S^{v} \right)}} = F_{0}},} & (56)\end{matrix}$where F₀ is a right-hand-side vector. In Eqn. (56) the saturations areexpressed in a vector notation. With the linearization of F(S) aroundS^(v), $\begin{matrix}{{{F\left( S^{v} \right)} + {\left( \frac{\partial F}{\partial S} \right)_{v} \cdot \left( {S^{v + 1} - S^{v}} \right)}} = {F_{0}.}} & (57)\end{matrix}$one obtains the iteration equation $\begin{matrix}{S^{v + 1} = {S^{v} + {\left( {F_{0} - F} \right)_{v} \cdot {\left( \frac{\partial F}{\partial S} \right)_{v}^{- 1}.}}}} & (58)\end{matrix}$

In addition, 0≦S_(j) ^(v+1)≦1 for all j is enforced after eachiteration.

In a 3-dimensional model of finite difference simulation, the gravityparameter, C_(g), is dependent on the flow direction, (i.e., C_(g)=0 forhorizontal flow F_(h), and C_(g)≠0 for vertical flow F_(v)). As aresult, two fractional flow curves have to be considered for horizontalF_(h) and vertical flow F_(v). The vertical flow F_(v) may become anon-monotonic curve for mC_(g)>1. The lighter fluid saturation (S₀)equation is solved and the other phase is computed by the saturationconstraint (S_(w)=1−S₀).

Consideration of Gravity

The effect of the additional consideration of gravity is potentiallyaltering the fractional flow curve F=f(S), Eqns. (31) and (32), so thatthose curves may contain a reflection point S*_(r). As described above,C_(g), Eqn. (30) is calculated and is used in Eqns. (24) and (25) inultimately establishing this fractional flow curve F(S).

In the event that there is no countercurrent flow, the fractional flowcurves will remain, as described in section I, without any reflectionpoint S*_(r). However, if there is countercurrent flow, then curves,such as seen in FIGS. 19A-B may contain a reflection point S*_(r).

Regions I and II may be defined as the regions bifurcated by thereflection point S*_(r) with Region I being the region below thereflection point, 0<S<S*_(r), and Region II being the region aboveS*_(r)<S<1.0.

In the event that a grid cell for which S^(v) is being evaluated is nextto a no flow boundary, the above-described method of section I isapplied without regard to any reflection point S*_(r). This is becausecountercurrent flow cannot occur under these conditions. Alternatively,if S^(v) in adjacent cells i and n are in different regions I and II,then the above-described method of section I is used to allow thesaturation S^(v+1) to move across the reflection point S*_(r).

In the event grid cell i is (1) not next to a no-flow boundary or else(2) cells i and n are in the same region, then S^(v+1) is determined forthe particular region in which S^(v) is located. For example, if S^(v)is in region I, then S^(v+1) will be maintained in region I. Iff″(S^(v))·f″(S^(v+1))≧0, then no underelaxation is needed to ensurestability. S^(v+1) will however be set to either 0 or S*_(r) in theevent S^(v+1) falls outside of region I. If f″(S^(v))·f″(S^(v+1))<0,i.e. on the opposite side of the inflection point S^(c) in that region,then underrelaxation will be used to insure stability. Again, if S^(v+1)falls outside of region I in which S^(v) is located, S^(v+1) will be setto one of the lower or upper boundary limits, 0 or S*_(r), prior to thecalculaton for the updated or underrelaxed value of S^(v+1) beingcalculated, i.e. 0.5(S^(v)+S^(v+1)).

Similarly, in region II, S^(v) and S^(v+1) are checked to see if S^(v)and S^(v+1) fall on the same side of the inflection point S^(c) inregion II. If so, no underrelaxation is needed for stability—only acheck is made to insure the S^(v+1) is within the range of Region II. IfS^(v) and S^(v+1) are on opposite sides of the inflection point S^(c) inthat region, i.e. f″(S^(v))·f″(S^(v+1))<0, then underrelaxation is used,i.e. S^(v+1)=0.5(S^(v)+S^(v+1)). Again, S^(v+1) must fall within regionII prior to doing the subsequent underrrelaxation calculation.

Capillary pressure is easy to incorporate into the above steps. {tildeover (C)}_(g) is simply set to C_(g)−C_(a), as will be described belowin Eqn. (65).

At this point S^(v+1) has been determined for this iteration v+1 forthis particular grid cell i. As described in Section I, a check forconvergence is performed on each of the grid cells of the model orsubmodel to see whether convergence has been achieved in all grid cellsso that a time step should be considered complete and no further Newtoniterations are necessary for the time step.

Algorithm for Updating S^(v)→S^(V+1) in Grid Cell

The following algorithm may be used for two-phase flow with gravity.

1. Check the directions of phase velocities and gravity force.

2. If the phase velocities are co-current, apply the regular stabilityanalysis;

-   -   Calculate the directional fractional flows and its second-order        derivatives (F″_(h) and F″_(v));    -   Calculate C_(g) for vertical flow (for horizontal flow,        C_(g)=0); and    -   Apply saturation limit 0≦S^(v+1)≦1.    -   If F″_(h)(S^(v+1))F″_(h)(S^(v))<0 or        F″_(v)(S^(v+1))F″_(v)(S^(v))<0, apply an under-relaxation for        saturation iteration: i.e., S^(v+1)=0.5(S^(v)+S^(v+1)).        3. If the phase velocities are counter-current;    -   Calculate the reflection saturation, S*_(r); and    -   apply saturation limit: if the cell is not adjacent to the        boundary or the neighboring cell doe not belong to the different        saturation domain, 0≦S^(v+1)≦S*_(r) or S*_(r)<S^(v+1)<1,else        0≦S^(v+1)≦1.    -   If F″_(h)(S^(v+1))F″_(h)(S^(v))<0 or        F″_(v)(S^(v+1))F″_(v)(S^(v))<0, apply an under-relaxation for        saturation iteration: i.e., S^(v+1)=0.5(S^(v)+S^(v+1)).        Capillary Pressure

In multiphase flow in porous media, fluids tend to have differentpressures due to the wettability to rock and the surface tension betweenfluids. This capillary pressure is often modeled by a simple function ofsaturation:p _(c)(S)=p ₂ −p ₁  (59)It is assumed that phase 1 is wetting and phase 2 is non-wetting. Ingeneral, the capillary pressure controls micro-scale flow in pores, butit plays a minor role in reservoir simulation scale. Therefore, it isnot expected that capillary pressure changes any characteristics of thestability calculation of transport equation.

Darcy's law may be written with capillary pressure: $\begin{matrix}{u_{1} = {{- \frac{{kk}_{r\quad 1}}{u_{1}}}\left( {{\nabla p_{1}} - {\rho_{1}g{\nabla D}}} \right)}} & (60) \\{u_{2} = {{- \frac{{kk}_{r\quad 2}}{u_{2}}}\left( {{\nabla p_{2}} - {\rho_{2}g{\nabla D}}} \right)}} & (61)\end{matrix}$

From Eqns. (59)-(61),u ₁ =f ₁(S)u−f ₃(S)(C _(g) −C _(a))  (62)u ₂ =f ₂(S)u+f ₃(S)(C _(g) −C _(a))  (63)where $\begin{matrix}{C_{a} = {\frac{k}{\mu_{1}}{\nabla p_{c}}}} & (64)\end{matrix}$

The transport equation does not change the characteristics due to thepresence of capillary pressure. The effect of capillary pressure can beincluded by modifying the gravity parameter C_(g). For the case of totalvelocity and gravity having different directions (i.e., phase 1 withC_(g)>0), capillary pressure reduces the gravity effect. The capillarypressure includes the spacial gradient of capillary pressure (∇p_(c)),and it is expected that this term stays much smaller than the viscousterm μ or gravity term C_(g). If there is any capillary force, the samestability algorithm may be used by modifying C_(g):{tilde over (C)} _(g) =C _(g) −C _(a).  (65)III. Algorithm for Three-Phase Flow with Gravity and Capillary Pressure

An unconditionally stable algorithm for three-phase flow with gravityand capillary pressure will now be described. In reservoir simulation,three-phase flow region is very limited in the model and more often thethree-phase region in the model occur is due to the combination of twotwo-phase flow. As a result, it is not clear that three phase relativepermeability measurements at laboratories can be directly applied forreservoir simulation. Furthermore, measurements of three-phase relativepermeability are technically very difficult, M. J. Oak, L. E. Barker,and D. C. Thomas, Three-phase relative permeability of Brea sandstone,J. Pet. Tech., pages 1057-1061, 1990. For given uncertainties inexperimental measurements and the scaling issues of grid-block relativepermeabilities in a reservoir model, a simple power law model ofrelative permeabilites for the stability analysis is applied. Thegeneral characteristics of fractional flow curves for three phase flowwill be examined. Based on the characteristics of the fractional flowcurves, an unconditionally stable algorithm has been designed forthree-phase flow with gravity and capillary pressure.

Three-Phase Flow

Three-phase flow in porous media is not well understood owing to thecomplexity of physical phenomena and practical difficulty in measuringthree-phase relative permeability. However, in reservoir simulation,three-phase flow in a grid cell is often seen that may occur mostly bythe combinations of two two-phase flows. Nevertheless, it entailsthree-phase relative permeabilities in reservoir simulation that may bean averaged value of two relative permeabilities for water-oil andoil-gas systems.

The modified Stone II method is often employed in many simulators inwhich oil relative permeability is modeled by a complicated functionalform: $\begin{matrix}{k_{ro} = {k_{rocw}\left\lbrack {{\left( {\frac{k_{row}}{k_{rocw}} + k_{rw}} \right)\left( {\frac{k_{rog}}{k_{rocw}} + k_{rg}} \right)} - k_{rw} - k_{rg}} \right\rbrack}} & (66)\end{matrix}$

Here, k_(r0), k_(rw), and k_(rg) are oil, water and gas relativepermeability, respectively. The k_(rocw) is oil relative permeability atconnate water saturation, k_(row) is oil relative permeability inoil-water system, and k_(rog) is oil relative permeability in oil-gassystem. The water and gas relative permeabilities are the same as thosegiven by a two phase system. The oil relative permeability of themodified Stone II method is found to become negative for some systems.Scant experimental evidence does not warrant employing a complexnonlinear functional form such as the Stone II method for three-phaserelative permeability. A. Heiba, Three-phase relative permeability,COFRC Technical Memorandum, 1987, proposed a straight-line interpolationof k_(row) and k_(rog). For given uncertainties in experimentalmeasurements and the scaling issues of grid-block relativepermeabilities in a reservoir model, a simple power law model for thestability analysis may be used. The analysis can be further extended forany reasonable three-phase relative permeability models (no negativevalues or non-monotonic curves).

Three-phase flow with gravity will first be examined. As in two-phaseflow stability analysis above, capillary pressure will later be includedby modifying a gravity number.

Darcy's law for three phase flow can be written as: $\begin{matrix}{u_{i} = {\frac{{kk}_{ri}}{u_{i}}{\left( {{\nabla p} - {\rho_{i}g{\nabla D}}} \right).}}} & (67)\end{matrix}$

Here, k is permeability, u is phase velocity, k_(r) is relativepermeability, μ is viscosity, D denotes the depth and ρ is density. Thesubscript i denotes phase properties (i.e., i=1, 2, 3 for water, oil andgas phases, respectively). From Eqn. (67), it can be readily derivedthat:u ₁ =f ₁(S ₁ ,S ₂ ,S ₃)u−f ₄(S ₁ ,S ₂ ,S ₃)C _(g1) −f ₅(S ₁ ,S ₂ ,S ₃)C_(g2)  (68)u ₂ =f ₂(S ₁ ,S ₂ ,S ₃)u+f ₄(S ₁ ,S ₂ ,S ₃)C _(g1) −f ₆(S ₁ ,S ₂ ,S ₃)C_(g3), and  (69)u ₃ =f ₃(S ₁ ,S ₂ ,S ₃)u+f ₅(S ₁ ,S ₂ ,S ₃)C _(g2) −f ₆(S ₁ ,S ₂ ,S ₃)C_(g3).  (70)

The total velocity is defined as:u=u ₁ +u ₂ +u ₃  (71)and the fractional parameters are given by $\begin{matrix}{f_{1} = \frac{k_{r\quad 1}}{k_{r\quad 1} + {m_{2}k_{r\quad 2}} + {m_{3}k_{r\quad 3}}}} & (72) \\{f_{2} = \frac{m_{2}k_{r\quad 2}}{k_{r\quad 1} + {m_{2}k_{r\quad 2}} + {m_{3}k_{r\quad 3}}}} & (73) \\{f_{3} = \frac{m_{3}k_{r\quad 3}}{k_{r\quad 1} + {m_{2}k_{r\quad 2}} + {m_{3}k_{r\quad 3}}}} & (74) \\{f_{4} = \frac{m_{2}k_{r\quad 1}k_{r\quad 2}}{k_{r\quad 1} + {m_{2}k_{r\quad 2}} + {m_{3}k_{r\quad 3}}}} & (75) \\{f_{5} = \frac{m_{3}k_{r\quad 1}k_{r\quad 3}}{k_{r\quad 1} + {m_{2}k_{r\quad 2}} + {m_{3}k_{r\quad 3}}}} & (76) \\{f_{6} = \frac{m_{2}m_{3}k_{r\quad 2}k_{r\quad 3}}{k_{r\quad 1} + {m_{2}k_{r\quad 2}} + {m_{3}k_{r\quad 3}}}} & (77)\end{matrix}$where m_(i) is the viscosity ratio with respect to the viscosity ofphase 1, μ₁/μ_(i). The velocity due to gravity force is characterized byvariables, C_(gi): $\begin{matrix}{{C_{g\quad 1} = {\frac{{kg}\left( {\rho_{1} - \rho_{2}} \right)}{\mu_{1}u_{c}}{\nabla D}}},} & (78) \\{{C_{g\quad 2} = {\frac{{kg}\left( {\rho_{1} - \rho_{3}} \right)}{\mu_{1}u_{c}}{\nabla D}}},{and}} & (79) \\{C_{g\quad 3} = {C_{g\quad 2} - C_{{g\quad 1},}}} & (80)\end{matrix}$where u_(c) is a characteristic velocity. Note that velocities u andu_(i) in (68)-(70) are nondimensionalized with respect to u_(c). Sincethere is a constraint for saturations:1=S ₁ +S ₂ +S ₃,  (81)only two equations in (68)-(70) are independent. Equations u₁ and u₂ arechosen for independent equations and also S₁ and S₂ are selected asindependent variables by substituting S₃ with 1−S₁−S₂.

It is instructive to examine fractional flow surfaces from Eqns. (68)and (69): $\begin{matrix}{F_{1} = \frac{u_{1}}{u}} & (82) \\{F_{2} = \frac{u_{2}}{u}} & (83)\end{matrix}$

An example is chosen with typical mobility ratios, m₂=0.2 and m₃=10 andquadratic relative permeabilities,k_(ri)=S_(i) ².  (84)

In FIGS. 20 and 21, the fractional flow surfaces are plotted for nogravity (C_(g1)=C_(g2)=0) and strong gravity (C_(g1)=5 and C_(g2)=10),respectively.

The fractional flow curves are monotonic for small C_(g1) and C_(g2),whereas it can become a multi-valued function for large C_(g1) andC_(g2). Note that depending on the directions of gravity force andvelocity, the fractional flow can be smaller than zero or larger thanone. This will cause numerical difficulty if the Newton nonlineariteration does not take into account the actual pattern of thesesurfaces.

From (68) and (69) the bounding curves can be derived where thefractional flow becomes zero:0=u−m ₂ C _(g1) k _(r2) −m ₃ C _(g2) k _(r3), and  (85)0=u+C _(g1) k _(r1) −m ₃ C _(g3) k _(r3.)  (86)

For quadratic relative permeabilities, Eqns. (85) and (86) become:$\begin{matrix}{{0 = {u - {m_{2}C_{g\quad 1}S_{2}^{2}} - {m_{3}{C_{g\quad 2}\left( {1 - S_{1} - S_{2}} \right)}^{2}}}},{and}} & (87) \\{0 = {u + {C_{g\quad 1}S_{1}^{2}} - {m_{3}{{C_{g\quad 3}\left( {1 - S_{1} - S_{2}} \right)}^{2}.}}}} & (88)\end{matrix}$Stability Analysis

As noted above, the stability issues of nonlinear transport problems liein the linearization of the flux function. The numerical schemes will beexamined as to whether they can converge from the initial saturationestimate, $\begin{matrix}{{{\lim\limits_{v\rightarrow\infty}{F\left( {S_{1}^{v},S_{2}^{v}} \right)}} = F^{*}},} & (89) \\{{{\begin{pmatrix}F_{1}^{v} \\F_{2}^{v}\end{pmatrix} + {\begin{pmatrix}\frac{\partial F_{1}}{\partial S_{1}} & \frac{\partial F_{1}}{\partial S_{2}} \\\frac{\partial F_{2}}{\partial S_{1}} & \frac{\partial F_{2}}{\partial S_{2}}\end{pmatrix}\begin{pmatrix}S_{1}^{v + 1} & {- S_{1}^{v}} \\S_{2}^{v + 1} & {- S_{2}^{v}}\end{pmatrix}}} = \begin{pmatrix}F_{1}^{*} \\F_{2}^{*}\end{pmatrix}},{and}} & (90) \\{{\begin{pmatrix}S_{1}^{v + 1} \\S_{2}^{v + 1}\end{pmatrix} = {\begin{pmatrix}S_{1}^{v} \\S_{2}^{v}\end{pmatrix} + {\begin{pmatrix}\frac{\partial F_{1}}{\partial S_{1}} & \frac{\partial F_{1}}{\partial S_{2}} \\\frac{\partial F_{2}}{\partial S_{1}} & \frac{\partial F_{2}}{\partial S_{2}}\end{pmatrix}^{- 1}\begin{pmatrix}F_{1}^{*} & {- F_{1}^{v}} \\F_{2}^{*} & {- F_{1}^{v}}\end{pmatrix}}}},} & (91)\end{matrix}$

In addition, 0≦S^(v+1)≦1 is enforced after each iteration.

Fractional Flow Surfaces

The fractional flow surfaces will now be examined to understand theirgeneral characteristics. In FIGS. 22, 23, and 24 the fractional flowisolines are plotted for three cases with no gravity (C_(g1)=C_(g2)=0),weak gravity (C_(g1)=0.5 and C_(g2)=1) and strong gravity (C_(g1)=5 andC_(g2)=10), respectively. Also applied is m₁=0.2 and m₃=10 for thisplot. If there is no gravity, the fractional flow for each phase movesin the same direction of the total velocity. There is no negative flowdomain in FIG. 22, where iso-fractional flow contours plotted from 0.1to 0.9 with the interval of 0.1. Note that a large area of small phaseone and two saturation has very weak fractional flows (<0.1). If thegravity is weak as in the second case, the domain of negative flow forphase 1 and 2 becomes smaller and the absolute magnitude of fractionalflow also becomes smaller. In the positive fractional flow in FIGS. 23and 24, the iso-fractional flow curves are contoured from 0 to 0.8 withthe interval of 0.2; but in the negative fractional flow domain,different contour values are used for each figure: (−0.2, −0.4, −0.6,−0.8, −1, −2) for FIG. 23A, (−0.05, −0.1 −0.15, −0.2, −0.25, −0.3) forFIG. 23B, (−0.02, −0.04, −0.06, −0.08, −1, −2) for FIG. 24A, and (−0.05,−0.1 −0.15) for FIG. 24B. Also plotted are the contours of u₁=0 (upperdotted line) and u₂=0 (lower dotted line). Due to the complex contourcurves in the negative flow domain, the Jacobian matrix of Eqn. (91) canbecome singular. The loci of singular points, where the determinant ofthe Jacobian matrix, Det, becomes zero, are also plotted with solidlines. A successful unconditionally stable algorithm should take intoaccount the singular points of Jacobian matrix because the generalNewton iteration cannot cross these singular lines.

The boundaries of zero fractional flow can be plotted for phase 1 (F₁=0)and phase 2 (F₂=0):L ₁(S ₁ ,S ₂):0=1−C _(g1) m ₂ S ₂ ² −C _(g2) m ₃(1−S ₁ −S ₂)², and  (92)L ₂(S ₁ ,S ₂):0=1+C _(g1) S ₁ ² −C _(g3) m ₃(1−S ₁ −S ₂)².  (93)

From Eqns. (93) and (94), the intersection point can be completed withS₂=0 axis. These points are denoted as S₁ ^(α1) from (92) and S₁ ^(α2)from (93), respectively. Two points are defined on these zero fluxcurves, (S₁ ^(β1)S₂ ^(β1)) and (S₁ ^(β2)S₂ ^(β2)), for better initialestimate of saturation that will be used in the algorithm describedlater.S ₁ ^(β1)=⅔S ₁ ^(α1) and S ₁ ^(β2)=⅓S ₁ ^(α2)  (94)L ₁(S ₁ ^(β1) ,S ₂ ^(β1)):0=1−C _(g1) m ₂(S ₂ ^(β1))² −C _(g2) m ₃(1−S ₁^(β1) −S ₂ ^(β1))², and  (95)L ₂(S ₁ ^(β2) ,S ₂ ^(β2)):0=1+C _(g1)(S ₁ ^(β2))² −C _(g3) m ₃(1−S ₁^(β2) −S ₂ ^(β2))².  (96)

A point is denoted in the middle of the domain for F*₁<0, F*₂<0, andDet>0, (S₁ ^(β3),S₂ ^(β3)) (i.e., (0.3, 0.3) in FIG. 24A).

Algorithm for Updating S^(v)→S^(v+1) in Three Phase Flow with Gravityand Capillary Pressure

-   1. Check the directions of total velocity and gravity force;-   2. Calculate C_(g1) and C_(g2);-   3. Calculate the boundaries of zero fractional flow for phase 1    (L₁(S₁,S₂)) and L₂(S₁,S₂);-   4. Calculate the initial points for saturation, (S₁ ^(β1),S₂ ^(β1))    and (S₁ ^(β2),S₂ ^(β2)) on L₁(S₁,S₂) and L₂(S₁,S₂), respectively;-   5. Check the sign of F*₁ and F*₂;-   6. From the initial estimate of saturations (S₁ ⁰ and S₂ ⁰),    calculate the fractional flow estimate (F₁ ⁰ and F₂ ⁰) and the    determinant Det of the Jacobian matrix;-   7. If F*₁≧0 and F*₂≧0;    -   (a) If the initial guess of the saturation does not yield the        positive condition, take the point at L₁(S₁ ^(β1),S₂ ^(β1)) as        the initial guess;    -   (b) Iterate the saturations based on Eqn. (91);    -   (c) If F₁ ^(v)<0 or F₂ ^(v)<0, set saturation at L₁(S₁ ^(β1),S₂        ^(β1));    -   (d) If S₁+S₂>1, normalize the saturation to ensure S_(i)≧0;    -   (e) If        ${{\frac{\partial^{2}F_{1}^{v + 1}}{\partial S_{1}^{2}}\frac{\partial^{2}F_{1}^{v}}{\partial S_{1}^{2}}} < 0},{or}$        $\frac{\partial^{2}F_{1}^{v + 1}}{\partial S_{2}^{2}}\frac{\partial^{2}F_{1}^{v}}{\partial S_{2}^{2}}$        <0, under-relax the saturation iteration, i.e.,        S^(v+1)=(S^(v)+S^(v+1))/2;-   8. If F*₁<0, F*₂≧0, and Det≧0, set the initial saturation estimate    at point at L₁(S₁ ^(β1),S₂ ^(β1)) as the initial guess;-   9. If F*₁<0, F*₂≧0, and Det<0, or F*₁<0, F*₂<0, and Det<0, set the    initial saturation estimate at point L₁(S₁ ^(β2),S₂ ^(β2)) as the    initial guess; and-   10. If F*₁<0, F*₂<0, and Det>0 set the initial saturation estimate    at point at (S₁ ^(β3),S₂ ^(β3)) as the initial guess.    Alternative Algorithm for Three Phase Stability

When three phase flow is modeled as a combination of two two-phase flowdomains (i.e., gas-oil and oil-water domains), the stability algorithmcan be modified as follows. The fraction of two-phase flow domain andtheir saturations (e.g., gas-oil and oil-water phases) are determinedfrom three phase relative permeability approximations (e.g., oilsaturation is uniform in the two domains or oil relative permeability isidentical in the two domains). The stability analysis for two-phase flowis applied for each domain in Newton-Raphson iterations. If the unstablesaturation change occurs in either one of the two domains, anunder-relaxation procedure is taken to warrant a stable convergence ofthe Newton-Raphson iterations.

Capillary Pressure

In multiphase flow in porous media, the fluids tend to have differentpressure due to their wettability to rock and interfacial tensionbetween fluids. This capillary pressure is often modeled by a simplefunction of saturation. In three-phase flow, we have two independentcapillary pressure curves:p _(c1)(S ₁ ,S ₂)=p ₂ −p ₁, and  (97)p _(c2)(S ₁ ,S ₂)=p ₃ −p ₁.  (98)The capillary pressure controls micro-scale flow in pores, but it playsa very minor role in reservoir simulation. It is not expected capillarypressure changes any characteristics of the stability calculation oftransport equation.

Darcy's law may be written with capillary pressure: $\begin{matrix}{{u_{1} = {{- \frac{{kk}_{r\quad 1}}{\mu_{1}}}\left( {{\nabla p_{1}} - {\rho_{1}g{\nabla D}}} \right)}},} & (99) \\{{u_{2} = {{- \frac{{kk}_{r\quad 2}}{\mu_{2}}}\left( {{\nabla p_{2}} - {\rho_{2}g{\nabla D}}} \right)}},{and}} & (100) \\{u_{3} = {{- \frac{{kk}_{r\quad 3}}{\mu_{3}}}{\left( {{\nabla p_{3}} - {\rho_{3}g{\nabla D}}} \right).}}} & (101)\end{matrix}$

Eqns. (97)-(101) yield:u ₁ =f ₁(S ₁ ,S ₂)u−f ₄(S ₁ ,S ₂)(C _(g1) −C _(a1))−f ₅(S ₁ ,S ₂)(C_(g2) −C _(a2))  (102)u ₂ =f ₂(S ₁ ,S ₂)u−f ₄(S ₁ ,S ₂)(C _(g1) −C _(a1))−f ₆(S ₁ ,S ₂)(C_(g3) −C _(a3))  (103)where $\begin{matrix}{{C_{a\quad 1} = {\frac{k}{\mu_{1}}{\nabla p_{c\quad 1}}}},} & (104) \\{{C_{a\quad 2} = {\frac{k}{\mu_{1}}{\nabla p_{c\quad 2}}}},} & (105) \\{C_{a\quad 3} = {C_{a\quad 2} - {C_{a\quad 1}.}}} & (106)\end{matrix}$

Note that the transport equation do not change the characteristics dueto the presence of capillary pressure. The effect of capillary pressurecan be included by the gravity parameter C_(g). For the case of totalvelocity and gravity have different directions (i.e., phase 1 withC_(g)>0), capillary pressure reduces the gravity effect. Since thecapillary pressure includes the spatial gradient of capillary pressure(∇p_(c)), this term is expectd to stay much smaller than the viscousterm u or gravity term C_(g). If there is capillary force, the samealgorithm of section II can be used by modifying C_(g):{tilde over (C)} _(gi) =C _(gi) −C _(a1).  (107)

1. A method for solving a non-linear S-shaped function F=f(S) which isrepresentative of a property S in a physical system, the S-shapedfunction having an inflection point S^(c) and a predetermined value of Fto which f(S) converges, the method comprising the steps of: a)selecting an initial value S^(v+1) as a potential solution to thefunction f(S); b) performing a Newton iteration (T) on the function f(S)at S^(v) to determine a next iterative value S^(v+1); c) determiningwhether S^(v+1) is located on the opposite side of the inflection pointS^(c) from S^(v); d) setting S^(v+1)=S^(l), a modified new estimate, ifS^(v+1) is located on the opposite side of the inflection point of f(S)from S^(v+1); e) determining whether S^(v+1) is within a predeterminedconvergence criteria; f) setting S^(v)=S^(v+1) if S^(v+1) is not withina predetermined convergence criteria; g) repeating steps (b)-(f) untilS^(v+1) is within the predetermined convergence criteria; and h)selecting S^(v+1) as a satisfactory solution S to the S-shaped functionF=f(S).
 2. The method of claim 1 wherein the modified new estimate S^(l)is set to S^(l)=αS^(v)+(1−α)S^(v+1), where 0<α<1.
 3. The method of claim2 wherein α=0.5.
 4. The method of claim 2 wherein 0.3<α<0.7.
 5. Themethod of claim 1 wherein the stability limit S^(l) is set to S^(c), theinflection point of the function f(S).
 6. The method of claim 1 whereinthe S-shaped function f(S) is a fractional flow function utilized in areservoir simulator.
 7. The method of claim 1 wherein the S-shapedfunction f(S) mathematically comports with the following mathematicalexpression:${f(S)} = \frac{\frac{S^{2}}{\mu_{1}}}{\frac{S^{2}}{\mu_{1}} + \frac{\left. \left( {1 - S} \right)^{2} \right)}{\mu_{2}}}$where S=saturation of a first phase in a two fluid flow problem; f(S)=afunction of saturation S; μ₁=viscosity of fluid in a first phase; andμ₂=viscosity of fluid in a second phase.
 8. The method of claim 1wherein S^(v+1) is within the predetermined convergence criteria ifε>|S^(v+1)−S^(v)| where ε is predetermined tolerance limit.
 9. Themethod of claim 1 wherein S^(v+1) is within the predeterminedconvergence criteria if ε>|F−f(S^(v+1))| where ε is predeterminedtolerance limit.
 10. The method of claim 1 further comprising settingS^(v+1) to an upper or lower bound S^(b) if S^(v+1) exceeds those boundsin step (b).
 11. The method of claim 1 wherein the S-shaped functionF=f(S) is characterized by one of a look-up table or a generalanalytical function.
 12. A program storage device readable by a machine,tangibly embodying a program of instructions executable by the machineto perform method steps for solving a non-linear S-shaped function,F=f(S), which is representative of a property S in a physical system,the S-shaped function having an inflection point S^(c) and apredetermined value of F to which f(S) converges, the method comprisingthe steps of: a) selecting an initial value S^(v) as a potentialsolution to the function f(S); b) performing a Newton iteration (T) onthe function f(S) at S^(v) to determine a next iterative value S^(v+1);c) determining whether S^(v+1) is located on the opposite side of theinflection point S^(c) from S^(v); d) setting S^(v+1)=S^(l), a modifiednew estimate, if S^(v+1) is located on the opposite side of theinflection point of f(S) from S^(v+1); e) determining whether S^(v+1) iswithin a predetermined convergence criteria; f) setting S^(v)=S^(v+1) ifS^(v+1) is not within a predetermined convergence criteria; g) repeatingsteps (b)-(f) until S^(v+1) is within the predetermined convergencecriteria; and h) selecting S^(v+1) as a satisfactory solution to theS-shaped function F=f(S).
 13. A method for updating a Newton's iterationon a fractional flow function f(S) associated with two-phase flow withgravity in a grid cell of a reservoir simulator, the method comprising:(a) checking the directions of phase velocities and gravity force in agrid cell; (i) if the phase velocities are co-current, applying aregular Newton stability analysis; (ii) calculating the directionalfractional flows F_(h) and F_(v) and its second-order derivatives(F″_(h) and F″_(v)); (iii) calculating C_(g) for vertical flow (forhorizontal flow, C_(g)=0); (b) applying a saturation limit 0≦S^(v+1)≦1;(i) if F″_(h)(S^(v+1))F″_(h)(S^(v))<0 or F″_(v)(S^(v+1))F″^(v)(S^(v))<0,apply an under-relaxation for saturation iteration; (c) if the phasevelocities are counter-current; calculating a reflection saturation,S*_(r); (i) applying a saturation limit: if the cell is not adjacent tothe boundary or the neighboring cell does not belong to the differentsaturation domain, 0≦S^(v+1)≦S*_(r) or S*_(r)<S^(v+1)<1, else0≦S^(v+1)≦1; and (ii) if F″_(h)(S^(v+1))F″_(h)(S^(v))<0 orF″_(v)(S^(v+1))F″^(v)(S^(v))<0, applying an under-relaxation forsaturation iteration.
 14. A method for updating a Newton's iteration ona fractional flow function f(S) associated with three-phase flow withgravity in a grid cell of a reservoir simulator, the method comprising:(a) check the directions of total velocity and gravity force; (b)Calculate C_(g1) and C_(g2); (c) calculainge the boundaries of zerofractional flow for phase 1; (L₁(S₁,S₂)) and L₂(S₁,S₂); (d) calculatingthe initial points for saturation, (S₁ ^(β1),S₂ ^(β1)) and (S₁ ^(β2),S₂^(β2)) on L₁(S₁,S₂) and L₂(S₁,S₂), respectively; (e) checking the signof F*₁ and F*₂; (f) From initial estimates of saturations (S₁ ⁰ and S₂⁰), calculate fractional flow estimate (F₁ ⁰ and F₂ ⁰) and thedeterminant of a Jacobian matrix; (g) if F*₁≧0 and F*₂≧0; (h) if theinitial guess of the saturation does not yield a positive condition,take the point at L₁(S₁ ^(β1),S₂ ^(β1)) as the initial guess; (i)Iterate the saturations based on Eqn. (91); (j) if F₁ ^(v)<0 or F₂^(v)<0, set saturation at L₁(S₁ ^(β1),S₂ ^(β1)); (k) if S₁+S₂>1,normalize the saturation to ensure S_(i)≧0; (l) if${{\frac{\partial^{2}F_{1}^{v + 1}}{\partial S_{1}^{2}}\frac{\partial^{2}F_{1}^{v}}{\partial S_{1}^{2}}} < 0},\quad{{or}\quad\frac{\partial^{2}F_{1}^{v + 1}}{\partial S_{2}^{2}}\frac{\partial^{2}F_{1}^{v}}{\partial S_{2}^{2}}}$<0, under-relax the saturation iteration; (m) if F*₁<0, F*₂≧0, and Det≧0set the initial saturation estimate at point at L₁(S₁ ^(β1),S₂ ^(β1)) asthe initial guess; (n) if F*₁<0. F*₂≧0, and Det<0, or F*₁<0, F*₂<0, andDet<0, set the initial saturation estimate at point L₁(S₁ ^(β2),S₂^(β2)) as the initial guess; and (o) if F*₁<0, F*₂<0, and Det>0 set theinitial saturation estimate at point (S₁ ^(β3),S₂ ^(β3)) as the initialguess.